Shocker—A Molecular Dynamics Protocol and Tool for Accelerating and Analyzing the Effects of Osmotic Shocks

The process of osmosis, a fundamental phenomenon in life, drives water through a semipermeable membrane in response to a solute concentration gradient across this membrane. In vitro, osmotic shocks are often used to drive shape changes in lipid vesicles, for instance, to study fission events in the context of artificial cells. While experimental techniques provide a macroscopic picture of large-scale membrane remodeling processes, molecular dynamics (MD) simulations are a powerful tool to study membrane deformations at the molecular level. However, simulating an osmotic shock is a time-consuming process due to slow water diffusion across the membrane, making it practically impossible to examine its effects in classic MD simulations. In this article, we present Shocker, a Python-based MD tool for simulating the effects of an osmotic shock by selecting and relocating water particles across a membrane over the course of several pumping cycles. Although this method is primarily aimed at efficiently simulating volume changes in vesicles, it can also handle membrane tubes and double bilayer systems. Additionally, Shocker is force field-independent and compatible with both coarse-grained and all-atom systems. We demonstrate that our tool is applicable to simulate both hypertonic and hypotonic osmotic shocks for a range of vesicular and bilamellar setups, including complex multicomponent systems containing membrane proteins or crowded internal solutions.


Figure S1 Shape analysis of steady-state simulations of different pumping rates.
A hypertonic shock was performed on POPC vesicles using pumping rates of 0.1% (slow) and 1% (fast) of the initial interior vesicle volume.After a 40% volume reduction, a 5 µs steady-state simulation was launched and the changes in shape for the two pumping rates are indicated.Each dot in the shape diagram represents the vesicle shape after 20 ns, corresponding to 10 evaluations of the vesicle volume.The color indicates the temporal evolution for the two pumping rates (white to dark red/blue).The snapshots of cut vesicles on the left and right show the final shape of the steady-state simulation for the fast and slow pumping rate, respectively.

Figure S2 Shape analysis of steady-state simulations for polymer-filled phase-separated vesicle.
A hypertonic shock was performed on a vesicle with a phase-separated DPPC/CHOL/DFPC membrane and filled with dextran and PEG.After a 35% volume reduction of the initial enclosed water volume, a 5 µs steady-state simulation was launched and the changes in shape are indicated.The color indicates the temporal evolution from white to dark green.The snapshot shows the equilibrated vesicle shape after 5 µs of steady-state simulation.

Figure S3
Oblicity analysis of steady-state simulations with and without transmembrane proteins.A hypertonic shock was performed on a POPC vesicle with and without aquaporin embedded into the membrane.After a 40% volume reduction, a 5 µs steady-state simulation was launched and the changes in oblicity S are calculated as described in the Supplementary Methods.A negative oblicity value indicates an oblate shape.The snapshots show the equilibrated vesicle shape after 5 µs of steady-state simulation from the side (cut) and top.

Figure S4 Histograms of water-water distances after relocation.
The center of mass distances of all relocated water particles to the closest other water particles in the bin for (a) a CG system and (b) an AA system.We applied pumping rates of 0.05% (yellow) and 1% (blue) for both systems, corresponding to 300 and 6000 water beads for the CG system or 30 and 600 water molecules for the AA system.The setup is a water box in which particles are moved from one side to the other.

Supporting Tables
Table S1 System composition and simulation details for all example systems.Shocker was used to perform hypertonic and hypotonic osmotic shock simulations for a variety of different systems.Simulation times varied between 200 ns and 2400 ns depending on the pumping rate.
a initial configuration b after pore closure

Supporting Methods
Analysis features of Shocker.The Shocker package has a few in-built analysis methods, which, if enabled, store useful properties of a vesicle each pumping cycle.In the first cycle, the inner and the outer membrane area as well as the vesicle volume are calculated.It is verified if the structure is whole or if it crossed the periodic boundaries of the simulation box.In the latter case, the smaller pieces are transferred to their original positions to make the structure whole again before further calculations.Subsequently, a point cloud is generated based on the lipid linker positions, in the case of Martini force field only the GL1 beads are considered.Alternatively, the user can provide the desired atom type for analysis.Subsequently, the Python library pyvista 1 was used to generate a triangulated surface of which the area and volume can be determined.This method allows for accurately calculating the area and volume of even the most irregularly shaped vesicles.In case solutes are present in the simulation box, it is possible to calculate the solute concentration inside and outside the vesicle.For this purpose, the vesicle volume from the previous step and the number of solute molecules are determined.The earlier calculated clusters are used to count the solutes in the inner and outer compartments.
For each pumping cycle, Shocker calculates the area of the inner leaflet (Ain), the area of the outer leaflet (Aex), and the interior volume (V) to calculate the relative volume v and the reduced area difference Δa (see Figure S6).Hereby, relative volume v is defined as and the reduced area difference Δa as where V0, Ain0 and Aex0 denote to the volume, area of the inner and outer leaflet at the beginning of the pumping simulation.A combination of these parameters quantifies the shape of a vesicle of the course of the pumping simulation. 2 Lastly, Shocker can determine the asphericity parameter Δ and oblicity parameter S.These values indicate the degree of asphericity and oblicity, respectively, of a set of points 3 and are calculated according to the eigenvalues of the gyration tensor with N being the total number of particles and riα is the α-th component of the position of particle i and α and β = x,y,z.The eigenvalues of T, denoted by λi, are the squares of the three principal radii of gyration.Thus,  .where  ̅ = 012 / is the average eigenvalue of the gyration tensor.The asphericity parameter Δ is 0 for a perfectly spherical set of points, while S is positive in the case of a prolate shape and negative for an oblate shape.

Figure
Figure S5 Evolution of potential energy during a pumping cycle.After 200 ps simulation time, water particles are moved for (a) a CG system and (b) an AA system as described in Figure S4.

Figure
Figure S6 Schematic representation of the different vesicle shape properties analyzed by Shocker.When a pumping simulation is initialized (left), Shocker determines the initial interior volume (V0), the area of the inner leaflet (Ain0) as well as the area of the outer leaflet (Aex0).At the end of each pumping cycle (right), the current volume (V), the area of the inner leaflet (Ain), and the area of the outer leaflet (Aex) are measured.With these values, Shocker can calculate the relative volume (v) and reduced area difference (Δa) which is used to describe shape changes of vesicles.